
library(dplyr)
library(RSwissMaps)
library(rjson)
library(ggeffects)
library(rstan)
library(ggplot2)
library(arm)
library(lme4)
library(rjson)
library(pxR)
library(performance)
library(lmeresampler)
library(broom.mixed)


##### Comparison with Masket model #####

make_data_party_canton <- function() {
  
  
  px_data <- read.px(file = url("https://dam-api.bfs.admin.ch/hub/api/dam/assets/32426125/master"), encoding = "UTF-8")
  data_muni <- as.data.frame(px_data)
  data_muni <- data_muni[data_muni$Ergebnis=="Ja in %",2:4]
  colnames(data_muni) <- c("DateName", "Place", "YesShare")
  
  data_json <- fromJSON(paste(readLines("https://www.pxweb.bfs.admin.ch/api/v1/fr/px-x-1703030000_101/px-x-1703030000_101.px"), collapse=""))
  muni_id <- as.data.frame(list(number_mun = data_json[["variables"]][[1]][["values"]],
                                Place = levels(data_muni$Place)))
  
  ## loads the ids for projects
  
  valid_items <- which(as.Date(substr(data_json[["variables"]][[2]][["valueTexts"]], 1, 10)) %in% as.Date(substr(levels(data_muni$DateName), 1, 10)))
  
  project_id <- as.data.frame(list(anr = data_json[["variables"]][[2]][["values"]][valid_items],
                                   DateName = levels(data_muni$DateName)))
  
  data_muni2 <- merge(data_muni, muni_id, by = "Place")
  data_muni3 <- merge(data_muni2, project_id, by = "DateName")
  
  data_muni3 <- data_muni3 %>% 
    mutate(place_type = ifelse(substr(Place, 1, 1)=="S", "Country",
                               ifelse(substr(Place, 1, 1)=="-", "Canton",
                                      ifelse(substr(Place, 1, 1)==">", "District",
                                             ifelse(substr(Place, 1, 1)==".", "Municipality", NA)))),
           id.place = ifelse(nchar(number_mun)<=4, number_mun,
                             ifelse(nchar(number_mun)>4, substr(number_mun, nchar(number_mun) - 3, nchar(number_mun)), NA)),
           year = as.numeric(substr(DateName, 1, 4)))
  
  data_muni3$anr <- as.numeric(data_muni3$anr)/10
  data_muni3$YesShare <- as.numeric(as.character(data_muni3$YesShare))
  
  
  ## Download the swissvotes data
  swissvotes <- read.csv(url("https://swissvotes.ch/page/dataset/swissvotes_dataset.csv"), sep = ";")
  
  ## Select columns and observation
  swissvotes <- swissvotes[swissvotes$anr>=193 & ### After 1965
                             swissvotes$anr<=664 & ###Before 2015 
                             swissvotes$rechtsform<=3, ## Only initivative, referenda mandatory and facultative
                           c(which(colnames(swissvotes)=="anr"),
                             which(colnames(swissvotes)=="datum"),
                             which(colnames(swissvotes)=="p.fdp"),
                             which(colnames(swissvotes)=="p.cvp"),
                             which(colnames(swissvotes)=="p.svp"),
                             which(colnames(swissvotes)=="p.gps"),
                             which(colnames(swissvotes)=="p.sps"),
                             which(colnames(swissvotes)=="rechtsform"),
                             which(colnames(swissvotes)=="volkja.proz"),
                             which(colnames(swissvotes)=="legislatur"))]
  
  valid_ids <- unique(swissvotes$anr)
  
  data_muni4 <- merge(data_muni3, swissvotes, by = "anr")
  
  data_canton <- data_muni4[data_muni4$place_type=="Canton",]
  data_mun <- data_muni4[data_muni4$place_type=="Municipality",]
  
  data_canton$i <- as.numeric(as.factor(data_canton$anr))
  data_canton$j <- as.numeric(as.factor(data_canton$number_mun)) 
  
  data_canton$CH_res <- rep(data_muni4[data_muni4$place_type=="Country",]$YesShare, each=26) ### Yes in percent at the national level
  
  data_canton$YesShare <- data_canton$YesShare/100 #### Divide the percent by 100 to have a 0 to 1 indicator
  data_canton$CH_res <- data_canton$CH_res/100
  
  data_canton$log_CH_res <- log(data_canton$CH_res/(1-data_canton$CH_res)) #### Compute the inverse logit of this value
  data_canton$log_YesShare <- log(data_canton$YesShare/(1-data_canton$YesShare)) 
  
  data_canton$y <- data_canton$log_YesShare - data_canton$log_CH_res
  
  data_canton <- data_canton %>% 
    mutate(date = as.Date(substr(DateName, 1, 10)),
           time = substr(date, 1, 4),
           t = as.numeric(as.factor(time)),
           t2 = as.numeric(as.factor(data_canton$legislatur)))
  
  data_mun$j <- as.numeric(as.factor(data_mun$number_mun)) 
  
  
  data_mun$i <- as.numeric(as.factor(data_mun$anr))
  
  data_mun$CH_res <- rep(data_muni4[data_muni4$place_type=="Country",]$YesShare, each=length(unique(data_mun$Place))) ### Yes in percent at the national level
  
  data_mun$YesShare <- data_mun$YesShare/100 #### Divide the percent by 100 to have a 0 to 1 indicator
  data_mun$CH_res <- data_mun$CH_res/100
  
  data_mun$log_CH_res <- log(data_mun$CH_res/(1-data_mun$CH_res)) #### Compute the inverse logit of this value
  data_mun$log_YesShare <- log(data_mun$YesShare/(1-data_mun$YesShare)) 
  
  data_mun$y <- data_mun$log_YesShare - data_mun$log_CH_res
  
  data_mun <- data_mun %>% 
    mutate(date = as.Date(substr(DateName, 1, 10)),
           time = substr(date, 1, 4),
           t = as.numeric(as.factor(time)),
           t2 = as.numeric(as.factor(data_mun$legislatur)))
  
  data_mun$j <- as.numeric(as.factor(data_mun$number_mun))
  
  
  for (i in 3:7) {
    if (i==3) {
      data_party <- swissvotes[,c(1, 2, 9)]
      data_party$pos <- swissvotes[,i]
      data_party$party <- substr(colnames(swissvotes)[i], 3, 5)
      data_party$year <- as.numeric(substr(swissvotes$datum, 7, 10))
      data_party$legislature <- swissvotes$legislatur
    } else {
      dta <- swissvotes[,c(1, 2, 9)]
      dta$pos <- swissvotes[,i]
      dta$party <- substr(colnames(swissvotes)[i], 3, 5)
      dta$year <- as.numeric(substr(swissvotes$datum, 7, 10))
      dta$legislature <- swissvotes$legislatur
      data_party <- rbind(data_party, dta)
    }
  }
  
  
  data_party2 <- as.data.frame(list(i = as.numeric(as.factor(data_party$anr)),
                                    j = as.numeric(as.factor(data_party$party)) + 26,
                                    t = as.numeric(as.factor(data_party$year)),
                                    t2 = as.numeric(as.factor(data_party$legislatur)),
                                    y = ifelse(data_party$pos==1, 1, ifelse(data_party$pos==2, 0, NA))))
  
  data_party3 <- as.data.frame(list(i = as.numeric(as.factor(data_party$anr)),
                                    j = as.numeric(as.factor(data_party$party)) + max(data_mun$j, na.rm = T),
                                    t = as.numeric(as.factor(data_party$year)),
                                    t2 = as.numeric(as.factor(data_party$legislatur)),
                                    y = ifelse(data_party$pos==1, 1, ifelse(data_party$pos==2, 0, NA))))
  
  
  data_canton2 <- as.data.frame(list(i = data_canton$i,
                                     j = data_canton$j,
                                     t = data_canton$t,
                                     t2 = data_canton$t2,
                                     y = ifelse(data_canton$YesShare<.5, 0, ifelse(data_canton$YesShare>.5, 1, NA))))
  
  
  
  data_mun2 <- as.data.frame(list(i = data_mun$i,
                                     j = data_mun$j,
                                     t = data_mun$t,
                                     t2 = data_mun$t2,
                                     y = ifelse(data_mun$YesShare<.5, 0, ifelse(data_mun$YesShare>.5, 1, NA))))
  
  
  
  
  data_party_canton <- rbind(data_party2, data_canton2)
  
  data_party_mun <- rbind(data_party3, data_mun2)
  
  
  return(list(data_party_canton=data_party_canton, 
              data_party_mun = data_party_mun))
  
}

data_party_canton_mun <- make_data_party_canton()

data_party_canton <- data_party_canton_mun$data_party_canton
data_party_mun <- data_party_canton_mun$data_party_mun

y <- data_party_canton$y
NAS <- which(is.na(y))
y <- y[-NAS]
i <- data_party_canton[-NAS,]$i
j <- data_party_canton[-NAS,]$j
t <- data_party_canton[-NAS,]$t
t2 <- data_party_canton[-NAS,]$t2
N <- length(y)
I <- max(i)
J <- max(j)
T <- max(t)
T2 <- max(t2)

data_IRT <- list(N=N, I=I, J=J, T=T, 
                 y=y, i=i, j=j, t=t)
data_IRT2 <- list(N=N, I=I, J=J, T=T2, 
                  y=y, i=i, j=j, t=t2)

y <- data_party_mun$y
NAS <- which(is.na(y))
y <- y[-NAS]
i <- data_party_mun[-NAS,]$i
j <- data_party_mun[-NAS,]$j
t <- data_party_mun[-NAS,]$t
t2 <- data_party_mun[-NAS,]$t2
N <- length(y)
I <- max(i)
J <- max(j)
T <- max(t)
T2 <- max(t2)


data_IRT_mun <- list(N=N, I=I, J=J, T=T, 
                 y=y, i=i, j=j, t=t)
data_IRT2_mun <- list(N=N, I=I, J=J, T=T2, 
                  y=y, i=i, j=j, t=t2)

model_party_canton <- stan(file = "Simple IRT.stan", data = data_IRT2, iter = 10000, warmup = 5000,
                           seed = 12345, cores=4, thin = 1)
save(model_party_canton, file = "model/Model IRT with parties and cantons.Rda")

model_party_mun <- stan(file = "Simple IRT.stan", data = data_IRT2_mun, iter = 1000, warmup = 500,
                           seed = 12345, cores=4, thin = 1)
save(model_party_mun, file = "model/Model IRT with parties and municipalities.Rda")

